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SUMMARY 


A  three-node  flat  shell  element  with  C®  rotation  fields  has  been  developed  for  analysis  of  arbitrary  composite 
shells.  The  element  may  consist  of  any  number  of  orthotropic  layers,  each  layer  having  different  material  properties 
and  angular  orientation.  The  formulation  includes  coupling  between  bending  and  extension,  which  is  essential  for 
analysis  of  unsymmetric  laminates.  Shearing  deflections  are  included,  since  laminated  and  sandwich  construction 
frequently  results  in  shear  stiffness  much  smaller  than  bending  stiffness.  Formulation  of  the  element  is  straight¬ 
forward,  and  calculation  of  its  stiffness  matrix  is  simple  and  fast.  Convergence  of  solutions  with  mesh  refinement 
is  uniform  for  both  thin  and  thick  shells  and  is  insensitive  to  element  shape,  although  not  as  rapid  as  some  other 
elements  that  lack  one  or  more  capabilities  of  the  newly  developed  element.  An  experimental  verification  of  the 
shell  element  is  reported  in  the  appendix. 

INTRODUCTION 

Significant  improvements  in  aircraft  flutter  speeds  and  other  performance  characteristics  are  possible  through  cre¬ 
ative  use  of  composite  materials,  as  has  been  pointed  out  by  Weisshaar  (1987).  Finite-element  analysis  of  structures 
of  this  type  requires  shell  elements  representing  layered  composites  that  may  be  unsymmetric  and  may  frequently 
use  honeycomb  or  foam  cores.  The  element  described  here  was  developed  for  this  purpose. 

Since  the  beginning  of  research  on  the  finite-element  method,  a  great  deal  of  effort  has  been  expended  on  devel¬ 
opment  of  shell  elements.  Many  papers  on  this  subject  have  been  published  (Gallagher,  1975;  1978).  However,  the 
perfect  shell  element  has  yet  to  be  invented. 

Triangular  flat-elements  having  displacements  and  rotations  at  the  comer  nodes  as  degrees  of  freedom  (dof)  are 
appealing  for  practical  reasons.  They  can  easily  model  arbitrary  shell  geometries  with  general  supports,  cutouts,  and 
beam  stiffeners.  These  elements  have  a  total  of  18  dof  (three  translations  and  three  rotations  at  each  node)  or  15  dof 
(three  translations  and  two  rotations),  depending  on  whether  the  rotation  about  the  normal,  which  has  zero  stiffness, 
is  included  as  a  degree-of-freedom. 

One  of  the  earliest  shell  elements  was  a  flat  triangular  element  developed  by  Melosh  (1966)  and  Martin  (1967). 
Versions  of  this  element  were  used  in  several  early  programs  including  ELAS,  SAMIS,  and  DYNAL.  At  that  time, 
it  was  the  only  element  that  used  the  transverse  shear  modulus,  and,  consequently,  it  was  the  only  element  suitable 
for  sandwich  plates  and  shells.  In  1971,  Utku  (1973)  developed  a  new  version  with  a  straightforward  treatment  of 
shearing  stiffness  which  is  theoretically  much  improved. 

The  element  developed  here  follows  the  derivation  by  Utku  (1973).  However,  integration  through  the  thickness 
is  done  stepwise,  required  for  representing  laminates,  and  two  additional  terms  of  the  energy  integrals  are  retained, 
since  they  do  not  vanish  for  unsymmetric  laminates. 

This  element  converges  at  about  the  same  rate  for  thick  or  thin  plates,  regardless  of  element  shape.  It  appears  to  be 
quite  robust  and  is  capable  of  correctly  modeling  behavior  of  general,  unsymmetric,  sandwich,  laminated  composite 
shells.  Convergence  is  slower  than  that  of  some  other  three-node  elements  that  do  not  include  shearing  deflections. 
Because  of  its  straightforward  formulation,  this  element  may  provide  a  basis  for  development  of  a  dynamic  element 
(Gupta,  1979),  which  could  result  in  a  significant  improvement  of  numerical  efficiency  in  analysis  of  problems  of 
aircraft  vibration  and  flutter. 

ELEMENT  FORMULATION 

The  triangular  element  is  developed  in  a  local  coordinate  system  with  origin  at  the  middle  surface  centroid,  and 
the  X  axis  is  parallel  to  the  element  1-2  edge,  as  shown  in  figure  1.  The  z  axis  is  perpendicular  to  the  plane  of  the 


element.  In  the  computer  implementation  of  the  element,  material  properties  for  a  layer  are  defined  with  respect  to 
a  material  coordinate  axis  which  may  be  oriented  at  an  angle  a  measured  from  the  element  x  axis  to  the  material 
1  axis.  The  material  3  axis  coincides  with  the  element  z  axis.  Displacements  in  the  x,  j/,  z  directions  are  denoted 
by  u,  V,  w.  The  usual  right-hand  vector  convention  is  used  for  directions  of  rotations,  which  are  denoted  by 
By,  and  6^.  The  rotations  and  By  are  rotations  of  lines  originally  perpendicular  to  the  middle  surface  of  the 
undeformed  element. 

Using  the  comma  convention  for  differentiation,  the  midsurface  strains  Co  are  given  by 


Co  = 


and  curvatures  are  given  by 
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In  this  derivation,  the  total  deflection  w  is  in  the  z  direction 


w  -  w'  +  wi^ 


where  w'  is  the  Kirchhoffian  part  of  the  deflection  and  w'*'  is  the  part  due  to  out-of-plane  shear. 
At  any  point  in  the  element,  the  strains  are  given  by 

€  =  €o  —  zX 


and  stresses  are  given  by 

a  =  De  =  Deo  —  DzX 

where  is  a  3  by  3  material  property  matrix  for  plane  stress  in  the  lamina. 
The  out-of-plane  shear  strains  are  given  by 
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and  the  corresponding  shear  stresses  are  given  by 


T  = 
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The  shear  moduli  in  D'  are  reduced  by  a  factor  of  5/6  to  account  for  nonuniform  distribution  of  shear  stress,  as  is 
done  in  many  other  programs. 

In  this  formulation,  a  straight  line,  which  is  normal  to  the  middle  surface  before  loading,  remains  straight,  but 
not  normal  to  the  deformed  middle  surface  after  loading.  The  slopes  of  the  middle  surface  are  and  w^y,  whereas 
the  slopes  of  the  line  that  was  originally  normal  to  the  middle  surface  become  —By  =  w^x  —  w^x  -  "^.y  ~  "^^y 

and  Bx  =  — 1',2  and  By  = 

The  total  strain  energy  in  the  element  is  given  by 


C7  =  ^  y  e^adAdz  2  J  I'^'^dAdz 
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The  energy  expression  can  be  expanded  into  five  terms,  as  follows: 


^  y  ^DeodAdz+  ^  J  DXz^  dAdz  +  ^  J  ')’^D'')dAdz+  ^  J  X’^ DcozdAdz  +  ^  J  ^DXzdAdz 


The  last  two  integrals  represent  coupling  between  middle-surface  deformations  and  bending  or  twisting  which  vanish 
if  the  laminate  is  symmetric  with  respect  to  the  middle  surface,  but  must  be  retained  in  this  development. 

In  integrating  through  the  thickness,  it  is  assumed  that  each  layer  is  a  different  homogeneous  orthotropic  material, 
but  none  of  the  matrices  depends  on  2.  For  example,  the  last  integral  in  the  aforementioned  energy  expression 
consequently  becomes 

^,DXdA 

^  t=i 

where  hi  is  the  z  coordinate  of  the  bottom  of  layer  i  and  n  is  the  number  of  layers. 

The  element  matrix  is  obtained  by  differentiation  of  the  energy  integral  with  respect  to  the  nodal  displacements, 
as  is  usually  done 

R- 

duiduj 

where  Ui,  uj  represent  all  of  the  nodal  displacements  and  rotations. 

To  proceed  with  development  of  the  element,  it  is  necessary  to  define  the  matrices  that  relate  strains  in  the  element 
to  nodal  displacements.  These  expressions  are  taken  from  the  paper  by  Utku  (1973),  and  the  details  of  the  derivation 
will  not  be  repeated  here.  The  nodal  displacements  in  the  x  direction  at  the  three  nodal  points  are  represented  by 


u  = 
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and  the  same  convention  is  used  to  represent  the  other  components  of  nodal  displacement. 

The  displacements  in  the  interior  of  the  element  u,v,6x,6y,  and  w*  are  each  interpolated  from  their  nodal  values 
by  the  same  linear  function.  For  example. 
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To  get  the  inteipolation  for  deflection  in  the  y  direction,  u  is  replaced  by  v  and  u  by  v,  and  so  forth.  In  this 
expression,  A  represents  area  of  the  element,  and  xi  and  yi  are  the  coordinates  of  node  1,  and  so  forth. 

It  follows  that  the  strains  on  the  middle  surface  Co  are  given  by 


eo=jj[MN] 
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and  the  curvatures  X  are  given  by 


where 


X=-[N-M] 
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and 


M  = 


y2  —  Vi  yj,  —  y\  y\  —  yi 
0  0  0 
13  —X2  XI  —  X3  X2  -  XI 


Development  of  the  expression  for  out-of-plane  shear  strain  7  in  terms  of  nodal  displacements  is  the  principal 
contribution  of  Utku’s  (1973)  paper.  It  is  given  by 

1  0  0  -2/1  0  0  XI  0  0  ■ 

0  10  0  -2/2  0  0  X2  0 

0  0  1  0  0  —2/3  0  0  X3 
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w 
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Perhaps  it  should  be  noted  that  some  terms  are  missing  from  the  matrices  which  are  the  end  result  in  Utku  (1973). 

In  the  computer  implementation  of  this  element,  the  matrices  relating  strains  to  nodal  displacements  are  calcu¬ 
lated  from  the  aforementioned  algebraic  expressions.  The  remaining  problems  are  calculated  numerically,  layer  by 
layer.  Since  the  strain-displacement  matrices  contain  only  constants,  integration  over  the  area  is  trivial.  Each  of 
the  five  energy  integrals  is  processed  to  yield  the  corresponding  submatrix  of  the  element,  stillhess  matrix,  and  the 
results  are  added  into  the  proper  location  in  the  element  stiffiiess  matrix.  The  program  then  goes  on  to  process  the 
next  layer  which  may  have  different  material  properties.  The  first  term  of  the  energy  integral,  for  example,  represents 
the  membrane  deformation  part  of  the  strain  energy  Um- 

Integrating  through  the  thickness,  this  becomes 
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where  hi  is  the  z  coordinate  of  the  bottom  of  layer  i  and  n  is  the  number  of  layers. 

Substituting  for  strains  in  terms  of  nodal  displacements  and  taking  second  partial  derivatives  of  the  energy  with 
respect  to  nodal  displacements,  results  in  the  submatrix  km,  which  represents  the  membrane  part  of  the  element 
stiffness  matrix.  ^ 

-  hi)lM  m'^DilM  N] 

Using  the  same  process 
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-  h^)[M  Nf  D[N  -M] 

i=l 

^^2  =  -  f^h[N  -MfDdMN] 

i=l 

Rotations  9x  and  9y  are  row  degrees  of  freedom  for  Displacements  u  amd  v  are  row  degrees  of  freedom  for 
kA2  •  The  element  stiffness  matrix  K  is 


■If  kfn  “t  JCff  +  +  ^j41  ^j42 

One  teim  in  each  of  the  aforementioned  summations,  which  represents  stiflhess  contribution  of  the  current  layer, 
is  computed.  Each  stiffness  coefficient  from  each  submatrix  is  added  into  the  appropriate  location  in  the  element 
stiffness  matrix.  The  program  then  goes  on  to  the  next  layer  until  all  layers  have  been  processed.  The  element 
stiffness  matrix  is  then  complete. 

ELEMENT  PERFORMANCE 

The  three-node  shell  element  developed  herein  should  prove  convenient  for  modeling  irregular  shapes  with 
cutouts,  as  well  as  for  analysis  of  unsymmetric  sandwich  composite  laminates  and  includes  out-of-plane  shear  de¬ 
formation.  In  comparison  with  other  elements,  it  may  be  tenned  LU71  (Laminated-Utku-1971).  Calculations  show 
that  the  LU71  gives  accurate  solutions  for  both  thick  and  thin  plates,  converges  uniformly  for  various  shapes,  and 
appears  robust.  Convergence  with  mesh  refinement  is  slower  than  that  of  the  two  elements  used  for  comparison,  but 
neither  of  those  elements  can  fill  the  need  for  which  the  LU71  was  developed. 

Figure  2  shows  one  of  the  test  problems  used  to  verify  the  element,  which  illustrates  the  coupling  of  bending 
and  in-plane  deflection.  Table  1  shows  calculated  deflections  for  four  cases,  all  with  the  same  loading.  These  results 
are  exact  and  can  be  easily  calculated  by  hand  because  strains  are  constant  throughout.  The  first  case  is  a  four-layer 
symmetric  laminate  with  fibers  only  in  the  0°  and  90°  directions.  The  second  case  represents  modeling  where  0°  and 
90°  elements  are  overlayed.  Results  are  identical.  The  third  case  is  an  unsymmetric  laminate.  Loads  applied  at  the 
middle  surface  cause  bending  and  out-of-plane  deflection.  The  fourth  case  is  identical  to  the  third,  except  that  the 
laminate  thickness  is  0.1  in.  instead  of  the  1.0  in.  used  in  the  first  three  cases.  Deflection  at  the  center  of  stiffness 
is  the  same  in  all  four  cases.  The  increase  in  the  value  of  the  deformation  in  the  X  direction  UX  for  the  third  and 
fourth  cases  is  due  to  rotation  and  the  calculation  of  the  UX  at  the  middle  surface. 

In  this  example.  Material  I  was  used  as  defined  by  Reddy  (1980)  with  Ei  =  10 .6  x  10  ^  psi. 

E1/E2  =  25 
Gi2/E2=0.5 
G2‘i/E2  =  0.2 
v\2  =  0 .25 

Figure  3  provides  a  convergence  comparison  of  the  LU71  element  with  the  DKT  element  (Batoz  and  others, 
1980),  which  is  type  53  in  ANSYS  (DeSalvo  and  Swanson,  1985),  and  with  Reddy’s  eight-node  isoparametric 
element  (Reddy,  1980)  for  a  thin  simply  supported  plate.  All  converge  to  the  same  result.  The  DKT  type  B  solution 
has  element  diagonals  at  —45°  rather  than  +45°,  and  convergence  is  slower,  but  stiff  very  good. 

The  convergence  in  figure  4  is  compared  with  that  in  figure  3  for  a  moderately  thick  plate.  The  DKT  element 
is  much  too  stiff  because  it  has  no  shearing  deflections.  The  LU71  element  converges  to  the  same  result  as  that  of 
Reddy’s  (1980)  isoparametric  element. 
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In  figures  3  and  4,  Reddy’s  result  is  plotted  as  being  equivalent  to  32  triangular  shell  elements,  since  it  is  estimated 
that  32  triangular  shell  elements  represent  about  the  same  computational  effort  as  the  four  isoparametric  quadrilateral 
elements  actually  used.  Reddy  uses  2  by  2  integration  in  his  examples,  but  for  some  8-node  elements,  this  can 
result  in  an  ill-conditioned  or  even  singular  matrix  in  structures  with  no  or  few  support  constaints.  In  application  to 
airframes,  it  is  believed  that  3  by  3  integration  would  be  required.  With  a  3  by  3  integration,  the  result  will  be  a  little 
less  accurate  than  that  shown  in  the  figures. 

The  DKT  element  (Batoz  and  others,  1980)  gives  outstanding  performance  in  most  bending  problems.  Its  con¬ 
vergence  rate  is,  however,  sensitive  to  boundary  conditions  and  element  shape.  For  membrane  deformations,  it  is 
a  constant  stress  triangle,  as  in  the  LU71.  Lack  of  shearing  deformations  is  its  greatest  weakness  in  analysis  of 
composites. 

Utku  (1973)  gives  seven  examples  comparing  the  performance  of  the  shell  element  formulation  used  here  with 
an  earlier  formulation  for  isotropic  homogeneous  plates.  The  examples  show  uniform  convergence  regardless  of 
element  shape.  The  formulation  used  here  converges  faster  than  the  earlier  one  for  equilateral  triangular  elements; 
but  the  convergence  is  slower  when  the  element  has  an  obtuse  angle. 

Table  2  shows  some  results  for  thin  plates  obtained  by  using  the  LU71.  The  example  is  a  48-in.^  homogeneous 
isotropic  plate  with  E  =  10 .6  x  10^  psi,  v  =  \/3  and  a  concentrated  load  of  16  lb  at  the  center.  AH  edges  were 
fixed.  One  quarter  of  the  plate  was  modeled  by  242  elements  in  a  regular  mesh.  Elements  were  right  triangles  with 
the  short  side  2.18  in.  long.  Since  the  LU71  includes  shear  deformation,  it  should  give  larger  deflections  for  the 
thicker  plates  but  the  same  deflection  for  very  thin  plates.  Table  2  shows  that  the  LU71  gives  results  as  expected  at 
the  plate  side  length/thickness  ratio  a/h  of  500.  Ata/h=  1000 ,  the  LU71  has  a  small  error,  for  thinner  plates,  the 
errors  become  greater. 

CONCLUDING  REMARKS 

The  three-node  shell  element  proves  to  be  convenient  for  modeling  irregular  shapes  with  cutouts,  and  is  also 
suitable  for  analysis  of  imsymmetric  sandwich  composite  laminates.  It  appears  to  be  robust  and  gives  accurate 
results  for  both  thick  and  thin  shells. 

Formulation  of  the  element  is  straightforward.  It  is  believed  that  this  formulation  can  be  used  as  a  basis  for  a 
dynamic  element  that  may  provide  much  improved  numerical  efficiency. 
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APPENDIX— COMPARISON  OF  EXPERIMENTAL  DEFLECTIONS  OF 
GRAPHITE  COMPOSITE  PLATES  WITH  FINITE  ELEMENT  CALCULATIONS 

Introduction 

Three  graphite-composite  plates  were  built  and  tested  to  verify  the  performance  of  the  LU71  finite  element.  The 
LU71  element  was  developed  during  1986  and  implemented  fortesting  in  a  program  called  LCSHEL  during  1987.  It 
is  a  three-node  C°  element  for  analysis  of  plates  and  shells  made  of  composite  materials  with  an  arbitrary  symmetric 
or  unsymmetric  stacking  sequence.  The  LU71  was  developed  to  meet  an  immediate  need  in  the  STARS  program 
(Gupta,  1984)  for  analysis  of  airframes  made  of  unsymmetric  laminated  composites  and  also  to  provide  a  foundation 
for  future  development  of  a  finite  dynamic  element  with  greatly  improved  numerical  efficiency.  The  experimental 
work  described  here  was  intended  to  verify  the  treatment  and  manipulation  of  material  properties  for  symmetric  and 
unsymmetric  laminates  and  the  overall  correcmess  of  analysis  by  using  the  element. 

Plate  Fabrication 

Three  graphite  fiber  composite  plates  were  fabricated.  The  graphite  fiber  was  supplied  in  a  12-in.-wide  roll 
of  unidirectional  fiber  with  a  few  cross-strands  of  glass  fiber  to  maintain  shape  and  fiber  spacing.  The  specified 
weight  of  this  material  was  0.033  Ib/ft^ ,  and  specified  tensile  strength  of  the  fiber  was  450,000  psi.  The  resin  used 
was  Hexcell  2410  with  2184  hardener.  Plates  were  placed  on  a  flat  surface,  vacuum  bagged,  and  cured  at  room 
temperature. 

Plate  A  was  a  four-layer  unidirectional  plate  9.25  by  12.0  in.  with  an  average  thickness  of  0.0303  in.  The  fiber 
fraction  was  0.572  by  weight.  Calculated  volume  fractions  were  Vf  =  0.474  fiber,  =  0.458  resin,  and  0.068 
voids.  This  plate  was  cut  into  strips  for  tension  tests  in  the  fiber  and  transverse  directions  for  determination  of 
material  properties. 

Plate  B  was  a  six-layer  unsymmetric  (0/0/0/90/90/90)  plate  12.0  by  12.0  in.  with  an  average  thickness  of 
0.04311  in.  The  fiberfraction  was  0.582  by  weight.  Calculated  volume  fractions  were  uy  =  0.494  fiber,  =  0.480 
resin,  and  0.0259  voids.  Thickness  was  measured  at  1-in.  grid  points  on  a  surface  plate  with  a  1/10,000-in.  dial  in¬ 
dicator.  The  average  thickness  was  0.04311  in.  with  a  standard  deviation  of  0.0058  in.  This  plate  was  trimmed  to 
11.3  by  11.3  in.  for  use  in  a  five-point  bend  test. 

Plate  C  was  an  eight-layer  symmetric  (0/0/90/90/90/90/0/0)  plate,  12.0  by  12.0  in.  with  an  average  thickness  of 
0.0549  in.  The  fiber  fraction  was  0.637  by  weight.  Calculated  volume  fractions  were  vy  =  0 .517  fiber,  =  0 .399 
resin,  and  0.084  voids.  Thickness  was  measured  at  1-in.  grid  points  on  a  surface  plate  with  a  1/10,000-in.  dial 
indicator.  The  average  thickness  was  0.0549  in.  with  a  standard  deviation  of  0.0046  in.  This  plate  was  trimmed  to 
1 1.3  by  11.3  in.  for  use  in  a  five-point  test. 

Material  Properties 

Material  properties  were  determined  from  tests  on  specimens  cut  from  plate  A.  Specimens  were  approximately 
1  in.  wide  by  9  in.  long.  Metal  plates  2.5  in.  long  were  glued  to  the  ends  of  the  specimens  (to  avoid  crushing),  leaving 
a  4-in.  test  length.  Biaxial  electric  resistance  strain  gages  were  glued  to  both  sides  of  the  specimens  and  all  gages 
were  read  individually. 

Material  properties  in  the  fiber  direction  were  consistent.  Unfortunately,  the  test  machine  for  measuring  trans¬ 
verse  properties  failed  to  work  properly,  and  some  data  were  lost.  Since  there  was  some  uncertainty  in  the  correct 
value  of  the  transverse  modulus,  the  material  properties  in  table  3  show  a  high  estimate  for  E2  (upper  bound)  and  a 
low  estimate  for  E2  Gower  bound). 
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Finite-element  calculations  were  performed  for  both  sets  of  properties  so  that  the  effect  of  this  uncertainty  could 
be  assessed.  Adjustments  in  properties  to  account  for  differences  in  fiber  volume  fractrons  were  done  using  relations 
from  Qiamis  (1984)  and  Jones  (1975). 

The  in-plane  shear  modulus  Gn  was  determined  from  a  four-point  bend  test,  supported  at  three  comers  and 
loaded  at  the  fourth  comer.  The  shear  modulus  is  determined  directly  as  given  by  Chandra  (1976).  When  limited  to 
the  linear  range  and  neglecting  self-weight,  this  simplifies  to 


in  which  o  is  plate  side  length,  h  is  plate  thickness,  P  is  load,  and  w  is  deflection  at  the  point  of  loading.  The 
value  of  Gi2  is  larger  than  would  be  expected  from  Chamis  (1984),  but  numerical  experiments  show  that  this  has  a 
negligible  influence  on  results.  This  larger  value  of  in-plane  shear  modulus  could  be  associated  with  slight  waviness 

in  the  fibers. 

Figure  5  shows  the  subscripting  convention  and  notation  for  elastic  constants. 

Plate  Bending  Tests 

Plates  B  and  C  were  tested  in  five-point  bending  with  supports  at  each  of  the  four  comers  and  a  concentrated  load 
at  the  center  (fig.  6).  Supports  were  set  0.35  in.  from  the  comers  of  the  plates  so  that  the  distance  between  supports 
was  10.6  in.  Loading  was  by  deadweights,  and  deflections  were  measured  by  a  microscope  with  a  resolution  of 
1/10,000  in. 

Plots  of  load  versus  deflection  data  are  shown  for  the  six-layer  unsymmetric  plate  in  figure  7  and  for  the  eight- 
layer  symmetric  plate  in  figure  8.  Lines  through  the  data  points  were  fitted  by  least  squares.  The  plots  show  good 
linearity  for  both  loading  and  unloading. 

Comparisons  to  Calculations 

The  finite-element  grid  used  for  calculations  is  shown  in  figure  9.  Calculations  were  done  with  the  HJ71  ele¬ 
ment  in  the  LCSHEL  program.  Additional  calculations,  not  included  here,  were  done  with  ANSYS,  and  agreement 
between  the  two  programs  was  good. 

Experimental  and  calculated  results  are  compared  in  table  4.  Since  there  was  some  uncertainty  about  the  values 
of  Et ,  calculations  were  done  separately  using  the  high  estimate  (upper  bound)  for  Ej  and  the  low  estimate  (lower 
bound)  for  Ez . 

Results  for  the  six-layer  unsymmetric  plate  (plate  B)  show  excellent  agreement.  Using  the  upper-bound  prop¬ 
erties,  the  calculated  result  is  2.3  percent  stiffer  than  the  experimental  result;  using  the  lower-bound  properties,  the 
calculated  result  is  2.2  percent  softer  than  the  experimental  result. 

Agreement  between  results  from  experiment  and  calculations  for  the  eight-layer  symmetric  plate  is  not  optimum, 
but  still  adequate.  Using  the  upper-bound  properties  results  in  a  calculation  5.7  percent  stiffer  than  the  experiment, 
while  using  the  lower-bound  properties  results  in  a  calculation  that  is  2.2  percent  stiffer  than  the  experiment.  The 
average  deviation  between  calculations  and  experiment  is  3.9  percent. 


Conclusion 

Comparison  of  results  from  plate  bending  experiments  and  finite-element  calculations  show  very  good  agree¬ 
ment,  and  support  the  conclusion  that  the  finite-element  program  is  using  composite  material  properties  correctly. 
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It  is  believed  that  the  resulting  discrepancies  are  due  primarily  to  imperfections  in  technique  for  fabricating  the 
composite  plates  and  performing  basic  property  measurements. 
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TABLE  1.  DEFLECTIONS  FOR  IN-PLANE  LOADING 


Stack 

sequence 

ux, 

10“'*  in. 

UY, 

10“^  in. 

uz, 

10  in. 

^S/i 

10“^  rad 

0/90/90/0 

-0.36207 

0.69629 

0 

0 

Superimposed 

0  and  90 

-0.36207 

0.69629 

0 

0 

0/90/90/0 

-1.0038 

1 .9303 

0.13898 

-Q. 21191 

0/0/90/90 

-10.038 

19.303 

13.898 

-27.797 

TABLE  2.  CENTER  DEFLECTION  OF  THIN  SQUARE  PLATES  WITH  FIXED  EDGES'* 


h,  in. 

a/h 

b/h 

Thin  plate 
theory,  in. 

LU71, 

in. 

Difference, 

percent 

0.96 

50 

2.27 

0.2180(10“^) 

0.2329  X 10“' 

-h6.8 

0.48 

100 

4.55 

0.1744(10“^) 

0.1840  xlO-2 

+5 .5 

0.096 

500 

22.7 

0.2180 

0.2231 

+2.3 

0.048 

1000 

45.5 

1.744 

1.684 

-3.4 

0.032 

1500 

68.2 

5.887 

5.287 

-10.0 

0.024 

2000 

90.9 

13.95 

11.56 

-17.2 

“Plate  side  length 

=  a  = 

48  in. 

Thickness  =  h 

Length  of  short  element  side  =  6  =  2 .18  in. 


TABLE  3.  MATERIAL  PROPERTIES 


Low  E2 
estimate 

High  El 
estimate 

Fiber  modulus,  Ef,l0^  psi 

30.13 

30.09 

Resin  modulus,  Em,  10^  psi 

0.2138 

0.2484 

0.3000 

0.3000 

0.1383 

0.0811 

Plate  B,  six  layers 

jBi,  10®  psi 

15.38 

15.38 

£2,10®  psi 

0.8463 

0.9784 

0.2181 

0.1886 

0.0120 

0.0120 

G'12, 10®  psi 

1.849 

1.849 

Plate  C,  eight  layers 

£;i,10®  psi 

15.64 

15.65 

f;2, 10®  psi 

0.7231 

0.8365 

1/12 

0.1911 

0.1868 

0.0088 

0.0100 

G'12, 10®  psi 

1.938 

1.938 

TABLE  4.  COMPARISON  OF  EXPERIMENTAL  MEASUREMENT  AND  FINITE-ELEMENT 
CALCULATIONS  OF  STIFFNESS  (LOAD/DEFLECTION  RATIO)  OF  GRAPHITE  COMPOSITE  PLATES 


Plate 

Properties 

Finite-element 
calculated  stiffness, 
Ib/in. 

Experimental 

stiffness, 

Ib/in. 

Difference, 

percent 

Plate  B 
(Six-layer) 

High  E2  estimate 

6.504 

6.359 

+2.3 

Low  Ez  estimate 

6.223 

6.359 

-2 .2 

Average 

6.364 

6.359 

0.0 

Plate  C 

(Eight-layer) 

High  Ez  estimate 

14.45 

13.67 

+5.7 

Low  Ez  estimate 

13.97 

13.67 

+2.2 

Average 

14.21 

13.67 

+3.9 

Figure  1.  Coordinate  systems. 


II 


stack  sequence 
[0/90/90/0] 
a  =  20  in. 
h  =  0.2  in. 
a/h  =  100 


O  Reddy’s  2  by  2  element 

-  LU71  program 

- ANSYS  DKT  element  type  B 

- ANSYS  DKT  element  type  A 


Stack  sequence 
[0/90/90/0] 
a  =  20  in. 
h  =  1  in. 
a/h  =  20 


O  Reddy’s  2  by  2  element 

-  LU71  program 

- ANSYS  DKT  element 


Normalized  15 
center 
deflection 


Symmetric 
Concentrated  load^ 


Symmetric 

5  ~  ^10  in.-] 

I  I  -J _ ^ ^ _ I _ ^ _ I 

0  25  50  75  100  125  150  175  200 

Number  of  elements  in  1/4  plate 

Figure  3.  Comparison  of  convergence  for  thin  plate. 


Normalized 

center 

deflection  ^2 


Symmetric 
Concentrated  load-\ 


Symmetric 
1-10  in.-| 


0  25  50  75  100  125  150  175  200 

Number  of  elements  in  1/4  plate 

Mb/ 

Figure  4.  Comparison  of  convergence  for  thick  plate. 


The  strains  for  plane  stress  due  to 
load  in  direction  1  are 


where  the  loading  direction  is 
denoted  by  left  superscript. 
Similarly,  the  strains  due  to  load 
in  direction  2  are 


therefore,  the  stiffness  matrix  for 
the  stress-strain  relations  for  plane 
stress  in  orthotropic  material  in 
terms  of  the  engineering  constants 
are 


El 

> 

> 

1 

V12E2 

1 

< 

ro 

_m 

1  -V12V21 

1 

^  1 

1 

< 

N> 

< 

E2 

1  -  V12  V21 

’  Qee  =  ®i2 

Stress  in  direction  1 
On 


0  £l 

0  £2 

Qee  ^12 


Stress  in  direction  2 


Qi2  0  £l  - 

02  =  0*12  ^22  ^ 

t.,2  0  0  Qee  j'i2_ 

Figure  5.  Material  property  definitions  (Jones,  1975,  pp.  39^6). 
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Figure  7.  Bending  test  of  six-layer  unsymmetric  plate.  Figure  8.  Bending  test  of  eight-layer  symmetric  plate. 


Figure  9.  Finite-element  mesh  used  for  calculation  of  stiffness  of  graphite  composite  plates. 
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